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Abstract 

We develop the equations of motion for full body models that describe the dynamics 
of rigid bodies, acting under their mutual gravity. The equations are derived using 
a variational approach where variations are defined on the Lie group of rigid body 
configurations. Both continuous equations of motion and variational integrators 
are developed in Lagrangian and Hamiltonian forms, and the reduction from the 
inertial frame to a relative coordinate system is also carried out. The Lie group 
variational integrators are shown to be symplectic, to preserve conserved quantities, 
and to guarantee exact evolution on the configuration space. One of these variational 
integrators is used to simulate the dynamics of two rigid dumbbell bodies. 



Key words: Variational integrators. Lie group method, full body problem 



Email addresses: tylee@umich.edu (Taeyoung Lee), mleok@umich.edu (Melvin Leok), 
nhm@engin.umich.edu (N. Harris McClamroch). 

^ This research has been supported in part by a grant from the Rackham Graduate School, 
University of Michigan. 

This research has been supported in part by NSF under grant ECS-0244977. 



Preprint submitted to Elsevier Science 



2 February 2008 



Nomenclat ure 



7j Linear momentum of the ith body in the inertial frame p. 12 

r Relative linear momentum p. 17 

Jj Standard moment of inertia matrix of the ith. body p. 9 

J(l^ Nonstandard moment of inertia matrix of the ith body p. 9 



Jr Standard moment of inertia matrix of the first body with respect to the p. 14 
second body fixed frame 

Jdj^ Nonstandard moment of inertia matrix of the first body with respect to p. 14 
the second body fixed frame 



nii Mass of the ith body p. 9 

m Reduced mass for two bodies of mass mi and m2, m = ^^^^^^ p. 17 

Mi Gravity gradient moment on the ith body p. 12 

fij Angular velocity of the ith body in its body fixed frame p. 9 

O, Angular velocity of the first body expressed in the second body fixed p. 13 
frame 

Ilj Angular momentum of the ith body in its body fixed frame p. 12 

n Angular momentum of the first body expressed in the second body fixed p. 17 
frame 

Ri Rotation matrix from the ith body fixed frame to the inertial frame p. 8 

R Relative attitude of the first body with respect to the second body p. 13 

R R = i?2) • ■ ■ ) -^n) P-9 

Vi Velocity of the mass center of the ith body in the inertial frame p. 12 

V Relative velocity of the first body with respect to the second body in p. 13 
the second body fixed frame 

V2 Velocity of the second body in the second body fixed frame p. 13 

Xi Position of the mass center of the ith body in the inertial frame p. 8 

X Relative position of the first body with respect to the second body ex- p. 13 
pressed in the second body fixed frame 

X2 Position of the second body in the second body fixed frame p. 13 

X X= (xi,X2,-- - ,Xn) P-9 
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1 Introduction 



1.1 Overview 



The full body problem studies the dynamics of rigid bodies interaeting under their mutual 
potential, and the mutual potential of distributed rigid bodies depends on both the position 
and the attitude of the bodies. Therefore, the translational and the rotational dynamics are 
coupled in the full body problem. The full body problem arises in numerous engineering and 
scientific fields. For example, in astrodynamics, the trajectory of a large spacecraft around 
the Earth is affected by the attitude of the spacecraft, and the dynamics of a binary asteroid 
pair is characterized by the non-spherical mass distributions of the bodies. In chemistry, 
the full rigid body model is used to study molccTilar dynamics. The importance of the full 
body problem is summarized in [7], along with a preliminary discussion from the point of 
view of geometric mechanics. 

The full two body problem was studied by Maciejewski [12], and he presented equations of 
motion in inertial and relative coordinates and discussed the existence of relative equilibria 
in the system. However, he does not derive the equations of motion, nor does he discuss the 
reconstruction equations that allow the recovery of the inertial dynamics from the relative 
dynamics. Scheeres derived a stability condition for the full two body problem [20], and 
he studied the planar stability of an ellipsoid-sphere model [21]. Recently, interest in the 
full body problem has increased, as it is estimated that up to 16% of near-earth asteroids 
are binaries [13]. Spacecraft motion about binary asteroids have been discussed using the 
restricted three body model [22] , [2] , and the four body model [23] . 

Conservation laws are important for studying the dynamics of the full body problem, 
because they describe qualitative characteristics of the system dynamics. The representa- 
tion used for the attitude of the bodies should be globally defined since the complicated 
dynamics of such systems would require frequent coordinate changes when using repre- 
sentations that are only defined locally. General numerical integration methods, such as 
the Runge-Kutta scheme, do not preserve first integrals nor the exact geometry of the full 
body dynamics [4]. A more careful analysis of computational methods used to study the 
full body problem is crucial. 

Variational integrators and Lie group methods provide a systematic method of constructing 
structure-preserving numerical integrators [4]. The idea of the variational approach is to 
discrctizc Hamilton's principle rather than the continuous equations of motion [17]. The 
numerical integrator obtained from the discrete Hamilton's principle exhibits excellent 
energy properties [3], conserves first integrals associated with symmetries by a discrete 
version of Noether's theorem, and preserves the symplectic structure. Many interesting 
differential equations, including full body dynamics, evolve on a Lie group. Lie group 
methods consist of numerical integrators that preserve the geometry of the configuration 
space by automatically remaining on the Lie group [5]. 

Moser and Vesselov [18], Wendlandt and Marsden [25] developed numerical integrators 
for a free rigid body by imposing an orthogonal constraint on the attitude variables and 
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by using unit quaternions, respectively. The idea of using the Lie group structure and 
the exponential map to numerically compute rigid body dynamics arises in the work of 
Simo, Tarnow, and Wong [24], and in the work by Krysl [8]. A Lie group approach is 
explicitly adopted by Lee, Leok, and McClamroch in the context of a variational integrator 
for rigid body attitude dynamics with a potential dependent on the attitude, namely the 
3D pendulum dynamics, in [9]. 

The motion of full rigid bodies depends essentially on the mutual gravitational potential, 
which in turn depends only on the relative positions and the relative attitudes of the bodies. 
Marsden et al. introduce discrete Euler-Poincare and Lie-Poisson equations in [14] and [15]. 
They reduce the discrete dynamics on a Lie group to the dynamics on the corresponding 
Lie algebra. Sanyal, Shen and McClamroch develop variational integrators for mechanical 
systems with configuration dependent inertia and they perform discrete Routh reduction 
in [19]. A more intrinsic development of discrete Routh reduction is given in [10] and [6]. 



1.2 Contributions 



The purpose of this paper is to provide a complete set of equations of motion for the full 
body problem. The equations of motion are categorized by three independent properties: 
continuous / discrete time, incrtial / relative coordinates, and Lagrangian / Hamiltonian 
forms. Therefore, a total of eight types of equations of motion for the full body problem 
are given in this paper. The relationships between these equations of motion are shown in 
Fig. 1, and are further summarized in Fig. 7. 



Continuous time (Sec. 3) 



Inertial coordinates (Sec. 3.1) 



Relative coordinates (Sec. 3.2) 



Lagrangian (3.1.1) 



Hamiltonian (3.1.2) 



Lagrangian (3.2.1) 



Hamiltonian (3.2.2) 



Discrete time (Sec. 4) 



Inertial coordinates (Sec. 4.1) 



Relative coordinates (Sec. 4.2) 



Lagrangian (4.1.1) Hamiltonian (4.1.2) Lagrangian (4.2.1) Hamiltonian (4.2.2) 



Fig. 1. Eight types of equations of motion for the full body problem 

Continuous equations of motion for the full body problem are given in [12] without any 
formal derivation of the equations. We show, in this paper, that the equations of motion for 
the full body problem can be derived from Hamilton's principle. The proper form for the 
variations of Lie group elements in the configuration space lead to a systematic derivation 
of the equations of motion. 
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This paper develops discrete variational equations of motion for the full body model follow- 
ing a similar variational approach but carried out within a discrete time framework. Since 
numerical integrators are derived from the discrete Hamilton's principle, they exhibit sym- 
plectic and momentum preserving properties, and good energy behavior. They are also 
constructed to conserve the Lie group geometry on the configuration space. Numerical 
simulation results for the interaction of two rigid dumbbell models are given. 

This paper is organized as follows. The basic idea of the variational integrator is given in 
section 2. The continuous equations of motion and variational integrators are developed in 

section 3 and 4. Numerical simulation results are given in section 5. An appendix contains 
several technical details that are utilized in the main development. 



2 Background 



2. 1 Hamilton's principle and variational integrators 



The procedure to derive the Euler-Lagrange equations and Hamilton's equations from 
Hamilton's principle is shown in Fig. 2. 



Configuration Space 










Lagrangian 


















Action Integral 
<S = JlfLiq,q)dt 










Variation 
S(S = £<S^ = 


Legendre transform. 
p = ¥L{qA) 








Euler-Lagrange Eqn. 

d dL dL r\ 
dt dq dq ~ ^ 


Hamilton's Eqn. 
q = Hp,p^ -Hq 



Configuration Space 
(%,%+i) xQ 



Discrete Lagrangian 

Ld{qk,<ik+i) 



Action Sum 



Variation 
= = 


Legendre transform. 

Pk = ^L{q,q)\k 








Dis. E-L Eqn. 


Dis. Hami 

Pk = -I 
Pk+i = I 


ton's Eqn. 



Fig. 2. Procedures to derive the continuous and discrete equations of motion 

When deriving the equations of motion, we first choose generalized coordinates q, and a 
corresponding configuration space Q. We then construct a Lagrangian from the kinetic and 

the potential energy. An action integral (5 = J^*-' L{q,q)dt is defined as the path integral 
of the Lagrangian along a time-parameterized trajectory. After taking the variation of the 
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action integral, and requiring it to be stationary, we obtain the Euler-Lagrange equations. 
If we use the Legendre transformation defined as 



p-5q = FL(g, q) ■ 5q 

- A 
de 



L{q,q + €6q), (1) 



e=Q 

where dq G TgQ, then we obtain Hamilton's equations in terms of momenta variables 
p = ¥L{q,q) G T*Q. These equations are equivalent to the Euler-Lagrange equations [16]. 

There are two issues that arise in applying this procedure to the full body problem. The 
first is that the configuration space for each rigid body is the semi-direct product, SE(3) = 
M'^(s)S0(3), where S0(3) is the Lie group of orthogonal matrices with determinant 1. 
Therefore, variations should be carefully chosen such that they respect the geometry of the 
configuration space. For example, a varied rotation matrix G SO (3) can be expressed as 

where e G M, and r] G so (3) is a variation represented by a skew-symmetric matrix. The 
variation of the rotation matrix SR is the part of i?^ that is linear in e: 

R^ = R + e5R + 0{e^). 

More precisely, 5R is given by 

d 



6R= , 

de 



R' = Rrj. (2) 





The second issue is that reduced variables can be used to obtain equations of motion 
expressed in relative coordinates. The variations of the reduced variables are constrained 
as they must arise from the variations of the unreduced variables while respecting the 
geometry of the configuration space. The proper variations of Lie group elements and 
reduced quantities are computed while deriving the continuous equations of motion. 



Generally, numerical integrators are obtained by approximating the continuous Euler- 
Lagrange equation using a finite difference rule such as = {qk+i — qk)/h, where qk 
denotes the value of q(t) at t = hk for an integration step size /i G M and an integer k. 
A variational integrator is derived by following a procedure shown in the right column of 
Fig. 2. When deriving a variational integrator, the velocity phase space (q, q) G TQ of the 
continuous Lagrangian is replaced by (qk,qk+i) & Q x Q, and the discrete Lagrangian 
is chosen such that it approximates a segment of the action integral 

Ld{qk,qk+i)~ / L{qk,k+iit),qk,k+i{t))dt, 
Jo 

where qk,k+iit) is the solution of the Euler-Lagrange equation satisfying boundary condi- 
tions gfe,jk+i(0) = qk and qk,k+i{h) = qk+i- Then, the discrete action sum 6^ = 1] Ld{qk, qk+i) 
approximates the action integral 0. Taking the variations of the action sum, we obtain the 
discrete Euler Lagrange equation 

Dg^^d(gfc-i, qk) + Bq^Ldiqk, qk+i) = 0, 
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where Yiq^Ld denotes the partial derivative of with respect to qk- This yields a discrete 

Lagrangian map F^^ : {qk-i,qk) '—^ ilkiQk+i)- Using a discrete analogue of the Legendre 
transformation, referred to as a discrete fiber derivative FLfj : Q x Q ^ T*Q, variational 
integrators can be expressed in Hamiltonian form as 

Pk = -'Dq^^Ldiqk, Qk+i), (3) 
Pk+i = ^qk+iLd{qk, Qk+i)- (4) 

This yields a discrete Hamiltonian map Fi^ : iqk,Pk) '-^ {qk+i,Pk+i)- A complete develop- 
ment of variational integrators can be found in [17]. 



2.2 Notation 



Variables in the inertial and the body fixed coordinates are denoted by lower-case and 
capital letters, respectively. A subscript i is used for variables related to the ith body. The 
relative variables have no subscript and the kth discrete variables have the second level 
subscript k. The letters x, v, Q and R are used to denote the position, velocity, angular 
velocity and rotation matrix, respectively. 



The trace of A G 



is denoted by 



tv[A] = 5^ [A], 



i=l 



where [A]ii is the i, ith element of A. It can be shown that 

tr[AS] = ti[BA] = tr[S^y4^] = ti[A^B'^] , 

3 

tr[A^i?] = [AUBU 

p,q=l 

tr[PQ] = 0, 



for matrices A,BE: M"^", a skew-symmetric matrix P 
matrix Q = Q'^ e M"^". 



(5) 
(6) 
(7) 

^ and a symmetric 



The map S -.R^ i-^ R^^^ is defined by the condition that S{x)y = x x y for x, y G R^. For 
X = {xi,X2,X3) G R^, S{x) is given by 



Six) 



-X3 X2 

xs —xi 

—X2 Xi 



It can be shown that 



S{xf = -S{x), 
S{xxy) = S{x)S{y) - S{y)S{x) = yx^ - xy^ , 



(8) 
(9) 
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S{Rx) = RS{x)R^, 
S{x)'^S{x) = [x^x) /3x3 — xx^ = tr[xx'^] 1^x3 — xx^, 

for x,yeR^ and R G S0(3). 

Using homogeneous coordinates, we can represent an element of SE(3) as follows: 

G SE(3) 



(10) 
(11) 



R X 
1 



for X G and R G S0(3). Then, an action on SE(3) is given by the usual matrix multi- 
plication in M*^^**. For example 



Ri xi 




i?2 X2 




R1R2 R1X2 + Xl 


1 




1 




1 



for xi,X2 G and Ri,R2 G S0(3). 

The action of an element of SE(3) on can be expressed using a matrix- vector product, 
once we identify with x {1} C M^. In particular, we see from 



R Xl 




X2 




RX2 + Xl 


1 




1 




1 



that X2 I— > Rx2 + Xl. 



3 Continuous time full body models 



In this section, the continuous time equations of motion for the full body problem are 
derived in inertial and relative coordinates, and are expressed in both Lagrangian and in 
Hamiltonian form. 

We define O — 616263 as an inertial frame, and O^^ — Ei^Ei^Ei^ as a body fixed frame for 
the zth body, Bi. The origin of the i\h body fixed frame is located at the center of mass of 
body Bi. The configuration space of the zth rigid body is SE(3) = (s) S0(3). We denote 
the position of the center of mass of Bi in the inertial frame by xi G M^, and we denote the 
attitude of Bi by Ri G SO (3), which is a rotation matrix from the ith body fixed frame to 
the inertial frame. 



3.1 Inertial coordinates 



Lagrangian: To derive the equations of motion, we first construct a Lagrangian for the full 
body problem. Given (x^, Ri) G SE(3), the inertial position of a mass element of Bi is given 
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by Xi + RiPi, where pi G denotes the position of the mass element in the body fixed 
frame. Then, the kinetic energy of Bi can be written as 

Using the fact that Jg Pidmi = and the kinematic equation Ri = RiS{^i), the kinetic 
energy Tj can be rewritten as 

Ti{xi,^li) = \ I \\xi\\^ + \\S{VLi)pi\\^ drrii, 

= Im \\xif + ^tv[s{ni)Jd,s{nif] , (12) 

where rrii G M is the total mass of Bi, Qi G M'^ is the angular velocity of Bj in the body fixed 
frame, and J^^ = Jg. pipjdrrii G M^^'^ is a nonstandard moment of inertia matrix. Using 
(11), it can be shown that the standard moment of inertia matrix Jj = S{pi)'^ S{pi)dmi G 
R^^3 jg related to J^. by 

Ji = tr[JrfJ hx3- Jd,, 

and that the following condition holds for any fij G 

SiJiQi) = S{ni)Jd, + JdiSi^i). (13) 

We first derive equations using the nonstandard moment of inertia matrix J^^, and then 
express them in terms of the standard moment of inertia Ji by using (13). 

The gravitational potential energy U : SE(3)"' i-* M is given by 

„ „ „ X 1 x-^ f f Gdmidrrij , 

C/(X1,X2,--- ,X„,i?i,i?2,--- ,i?n =-7^ V / / II ^ „ ^, 14 

2 JBi JBj + RiPi - Xj - RjPjW 

where G is the universal gravitational constant. 

Then, the Lagrangian for n full bodies can be written as 

n 

L(x,x,R,n) = ll^if + ^^AS{^i)JdiS{nif] - C/(x,R), (15) 

i=l 

where bold type letters denote ordered n-tuples of variables. For example, x G (M^)", 
R G S0(3)", and n G (M^)" are defined as x = (xi,X2,-- - R = (i?i,i?2,--- ,Rn), 

and Q, = {0,1,^2, ■ ■ ■, ^n), respectively. 

Variations of variables: Since the configuration space is SE(3)", the variations should be 
carefully chosen such that they respect the geometry of the configuration space. The vari- 
ations of Xi : [to,tf] and Xi : [to,tf] i— are trivial, namely 

xj = Xi + e6xi + 0{e^), 
xl = Xi + eSxi + O(e^), 
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where 5xi : [to,tf] i— > R^, Sxi : [to,tf] i— > vanish at the initial time to and at the final 
time tf. The variation of Ri : [to,tf] S0(3), as given in (2), is 



SRi = RiTji, 



(16) 



where rji : [to,tf] i— 0o(3) is a variation with values represented by a skew-symmetric 
matrix (rjf = —rji) vanishing at to a^nd tj. The variation of can be computed from the 
kinematic equation Ri = RiS{0,i) and (16) to be 



s(«n.) = - 



Rf^ R\ = SRJ Ri + rJSRi, 



(17) 



3.1.1 Equations of motion: Lagrangian form 

If we take variations of the Lagrangian using (16) and (17), we obtain the equations of 
motion from Hamilton's principle. We first take the variation of the kinetic energy of Bi. 



6Ti 



de 



Ti(xi + eSxi, Qi + e6Qi), 



e=Q 



= mixJSxi + -tT[S{Sni)Jd^S{nif + S{ni)Jd,SiSnif] . 
Substituting (17) into the above equation and using (5), we obtain 
STi = mixjdxi + ^tr[-?7i {J<i.6'(rii) + 5(0^) JdJ] 



+ -tirii {Si^^i) {Jd,s{ni) + si^i)Jd,} - UdM^i) + s{ni)Jd,}si^}i)}] 



Using (9) and (13), 6Ti is given by 



1 



STi = miXi Sxi + -ti[- fiiS{Ji9.i) + r]iS{9.i x JjJ^j)] . 



(18) 



The variation of the potential energy is given by 

?7(x + e(5x,R + e(5R), 



e=0 



where 5x = {6x1,6x2, • • • , 6xn), 6R = {5Ri, 6R2, • ■ ■ , 6Rn). Then, 6U can be written as 



n / 3 



dU 



dU 



(19) 



where [A\pq denotes the p, gth element of a matrix A, and ^ are given by [^\p 
a^^' [^^ip? ~ d'^] ' respectively. The variation of the Lagrangian has the form 



5L = Y^ 5Ti - 6U, 



(20) 



i=l 
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which can be written more exphcitly by using (18) and (19). 
The action integral is defined to be 

/•*/ 

&= L(x,x,R,n)dt 

J to 

The variation of the action integral can be written as 

" ^ rtf du^ if r 

50 = ^ y mixfSxi - ^ Sxi + -tr -'qiS{Ji9.i) + r/i |S'(fJj X JjO^) + 2RI 
Using integration by parts, 

n 

6e = Y^ 



dU 



(21) 



dt. 



niiX^ 5xi 



i=l 



tf 1 
to 



*^ /■*/ 1 

+ / -rriiXi 6xi + -tr 
to •'to ^ 



r]iS{JiQ.i 



dt 



+ 2^ / -— Sxi + -ti 



=1 •'*o 



Using the fact that Sxi and r/j vanish at to and t/, the first two terms of the above equation 
vanish. Then, S(S is given by 



" rt 



i=l •^*o 



. T ) - dU] 1 
-dXi i rriiXi + ^ | + 



r/i <^ S'( Jifii + X Ji^lj) + 2i?; 



dt. 



From Hamilton's principle, (5(5 should be zero for all possible variations 5xi : [tQ^tf] 1— >■ 
and rji : [to;^/] ^ so (3). Therefore, the expression in the first brace should be zero. 
Furthermore, since rji is skew symmetric, we have by (7), that the expression in the second 
brace should be symmetric. Then, we obtain the continuous equations of motion as 



f^iXi — 



du_ 

dxi' 



dU 
dRi 



gTjT 

S{Ji(li + niX JiQif + 2— R. 



dRi 



(22) 



Using (8), we can simplify (22) to be 



dRi • dRi 

Note that the right hand side expression in the above equation is also skew symmetric. 
Then, the moment due to the gravitational potential on the ith body, Mj € is given 
by S{Mi) = Ri — Rf-§^- This moment can be expressed more explicitly as the 
following computation shows. 



S{Mi) 



dU 



Ri — Ri 



dU 
dRi' 



T T T 
"^rii '^ri2 '^ri-j, 



T T T 
h 12 is 



Uri 



Ur 



Tt3 
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where ri^,Urip G M^^^ are the pth row vectors of Ri and respectively. Substituting 
X = rj^, y = ujj^^ into (9), we obtain 

S{Mi) = Sin^ X Urii) + 5(ri2 X 14^2) + S(ri^ x 1*^3), 

= 'S'(rjj X 'u^jj + X u^jj + ''is x Uj-i^), (23) 

Then, the gravitational moment Mj is given by 

Mi = Til X Urii + X tiria + fia X ^^ris- (24) 



In summary, the continuous equations of motion for the full body problem, in Lagrangian 
form, can be written for z G (1, 2, • • • , n) as 

1 dU 

Vi = 25 

Jitli + Hi X Ji^i = Mi, (26) 
Xi = Vi, (27) 
Ri = RiSi^i), (28) 

where the translational velocity Vi G M? is defined as Vi = Xi. 



3.1.2 Equations of motion: Hamiltonian form 

We denote the linear and angular momentum of the ith body Bi by 7^ G M^, and Xlj G M^, 
respectively. They are related to the linear and angular velocities by 7^ = miVi, and Ilj = 
JjOj. Then, the equations of motion can be rewritten in terms of the momenta variables. 
The continuous equations of motion for the full body problem, in Hamiltonian form, can 
be written for i G (1, 2, • • • , n) as 

Ui + n^xUi = Mi, (30) 

Xi = —, (31) 
mi 

Ri = RiSiQi). (32) 



3.2 Relative coordinates 



The motion of the full rigid bodies depends only on the relative positions and the relative 
attitudes of the bodies. This is a consequence of the property that the gravitational po- 
tential can be expressed using only these relative variables. Physically, this is related to 
the fact that the total linear momentum and the total angular momentum about the mass 
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center of the bodies are conserved. Mathematically, the Lagrangian is invariant under a 
left action of an element of SE(3). So, it is natural to express the equations of motion in 
one of the body fixed frames. In this section, the equations of motion for the full two body 
problem are derived in relative coordinates. This result can be readily generalized to the 
n body problem. 

Reduction of variables: In [12], the reduction is carried out in stages, by first reducing 
position variables in M^, and then reducing attitude variables in S0(3). This is equivalent 
to directly reducing the position and the attitude variables in SE(3) in a single step, 
which is a result that can be explained by the general theory of Lagrangian reduction 
by stages [1]. The reduced position and the reduced attitude variables are the relative 
position and the relative attitude of the first body with respect to the second body. In 
other words, the variables are reduced by applying the inverse of {R2,X2) G SE(3) given 
by {R2, —R2X2) £ SE(3), to the following homogeneous transformations: 



rJ 


-R2X2 


( 


Ri xi 




R2 X2 


) 





1 




1 




1 






i?Ji?2 R2ix2 - X2) 

1 



(33) 



This motivates the definition of the reduced variables as 



X ■ 
R 



R2{xi - X2), 

R2R1, 



(34) 
(35) 



where X G is the relative position of the first body with respect to the second body 

expressed in the second body fixed frame, and R G SO (3) is the relative attitude of the 
first body with respect to the second body. The corresponding linear and angular velocities 
are also defined as 



V = R^ixi-X2), 

f2 = i?r2i. 



(36) 
(37) 



where 1/ G represents the relative velocity of the first body with respect to the second 
body in the second body fixed frame, and G M"^ is the angular velocity of the first 
body expressed in the second body fixed frame. Here, the capital letters denote variables 
expressed in the second body fixed frame. 

For convenience, we denote the inertial position and the inertial velocity of the second 
body, expressed in the second body fixed frame by X2, F2 G M^: 



X2 
V2 



R2X2, 
R2X2. 



(38) 
(39) 



Reduced Lagrangian: The equations of motion in relative coordinates are derived in the 
same way used to derive the equations in the inertial frame. We first construct a reduced 
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Lagrangian. The reduced Lagrangian I is obtained by expressing the original Lagrangian 
(15) in terms of the reduced variables. The kinetic energy is given by 



= ^m,\\v + V2f + ^m2\\V2f + ^tT[s{n)Jd^sinf] + ^tT[s{n2)^^^^ , 

where (10) is used, and J^^ G R-^^-^ is defined as J^j^ = RJ^^R^, which is an expression 
of the nonstandard moment of inertia matrix of the first body with respect to the second 
body fixed frame. Note that J^^ is not a constant matrix. Using (10), it can be shown that 
also satisfies a property similar to (13), namely 

s{jRn) = smja^ + Jd^sin), (40) 

where Jr = RJ\R?^ G R^^^ is the standard moment of inertia matrix of the first body with 
respect to the second body fixed frame. 

Using (14), the gravitational potential U of the full two bodies is given by 

U{xl,X2,Rl,R2) = - / T] ^ [7, 

and it is invariant under an action of an element of SE(3). Therefore, the gravitational 
potential can be written as a function of the relative variables only. By applying the inverse 
of {R2,X2) G SE(3) as given in (33), we obtain 

U{xi,X2,Ri,R2) = U{R^{xi-X2),0,R^Ri,l3y,3), 

Gdmidm2 



Bi Jb2 W^Ii^i - X2) + R2R1P1 - ^3x3P2|| ' 
Gdmidm2 



iBi Jb2 + -^Pi ~ P2II 
= U{X, R). 

Here we abuse notation slightly by using the same letter U to denote the gravitational 
potential function of the relative variables. 

Then, the reduced Lagrangian I is given by 



liR,X,^,V,^2,V2) = -mi\\V + V2f + -m2\\V-2''^ 



2 " " 2 

+ ^tv[s{n)Jd^sinf] + ltis{n2)Jd2S{n2f] - u{x, r). (4i) 



Variations of reduced variables: The variations of the reduced variables must be restricted 
to those that can arise from the variations of the original variables. For example, the 
variation of the relative attitude R is given by 

R^Rl, 

e=0 
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Substituting (16) into the above equation, 

SR = -r/2i?2 i?i + RlRim, 
= — 772-R + r]R, 

where a reduced variation 77 G 50 (3) is defined as ry = RrjiR^ . The variations of other 
reduced variables can be obtained in a similar way. The detailed derivations are given in 
A.l, and we summarize the results as follows: 

SR = r]R- r]2R, (42) 

dx = x- mx, (43) 

S{Sn) =rt- S{n)rt + rtS{n) + S{n)r)2 - V2S{n) + s{n2)v - VS{^2), (44) 
SV = x + S{n2)x - V2V, (45) 

s{sn2) = rj2 + s{n2)7]2-V2S{n2), (46) 

SV2 = X2 + S{n2)x2 - r]2V2, (47) 

where X:X2 '■ [to,tf] i— >■ and 77,7/2 : [io,i/] '-^ so(3) are variations that vanish at the end 
points. These Lie group variations are the key elements required to obtain the equations 
of motion in relative coordinates. 



3.2.1 Equations of motion: Lagrangian form 

The reduced equations of motion can be computed from the reduced Lagrangian using 
the reduced Hamilton's principle. By taking the variation of the reduced Lagrangian (41) 
using the constrained variations given by (42) through (47), we can obtain the equations 
of motion in the relative coordinates. 

Following a similar process to the derivation of STi,SU as in (18) and (19), the variation 
of the reduced Lagrangian SI can be obtained as 

SI = [miiV + V2)] - x^ [^1^2 x{V + V2)] 

+ [miiy + V2) + 7712^2] - xl [min2 X (y + V2) + m29.2 X V2] 



+ hi[-'qS{jR^l) + 775(02 X JrQ)] + ^tr[-7725( J2O2) + 772^(^2 x J2O2) 



- X 



dX 









+ tr 




+ tr 



(48) 



where we used the identities (9), (13) and (40), and the constrained variations (42) through 
(47). 



The action integral in terms of the reduced Lagrangian is 



Jto 



i{R,x,n,v,n2,V2)dt. 



(49) 
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Using integration by parts together with the fact that x^XitV and 772 vanish at Iq and tf, 
the variation of the action integral can be expressed from (48) as 



to 



5& = - <rni{V + V2)+min2X {V + V2) + 



dU 
dX 



dt 



dt 



+ 2 



'to 

1 r'f 



tr 



to 



7] { s{{jRn) + n2x Jnn) - 2r— 



dt 



+ 2 



1 rf 



tr 



to 



V2 <! S{J2n2 + X 72^2) + 2X— + 2R— 



dt. 



From the reduced Hamilton's principle, 6<5 = for all possible variations X)X2 '■ [to,tf] 

and 77,772 : [io;*/] '-^ 50(3) that vanish at to and tf. Therefore, in the above equation, 
the expressions in the first two braces should be zero and the expressions in the last two 
braces should be symmetric since 77, 7^2 are skew symmetric. Then, we obtain the following 
equations of motion, 

mi{V + V2) + min2x{V + V2) = - — , (50) 

m2V2 + 771,2^^2 ^ ^2 = (51) 

s{{jRn) + n2x JrU) = -s{M), 

8U BU 

S{J2tl2 + ^2>^J2^2)=Q^X^-X— +S{M), (52) 

where M G is defined by the relation S{M) = — R^^ ■ By a procedure analogous 

to the derivation of (23), M can be written as 

M = ri X Ur-^ + r2 X u^j + r^x u^^ , (53) 

where rp,Urp G are the pth column vectors of R and respectively. 

Equation (50) can be simplified using (51) as 

mi + 777-2 dU 



V + n2xv = — 



77117772 dX 



For reconstruction of the motion of the second body, it is natural to express the motion of 
the second body in the inertial frame. Since V2 = Rlx2 + Rlx2 = -S{^l2)V + Rlv2, (51) 
can be written as 

Equation (52) can be simplified using the property ^X"^ — X^^ = S{X x ^) from (9). 
The kinematics equations for R and X can be derived in a similar way. 
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In summary, the continuous equations of relative motion for the full two body problem, in 
Lagrangian form, can be written as 



(Jr^) + ^2X Jr^ = -M, 
dU 
oX 

X + ^2XX = V, 
R = S{n)R- S{n2)R, 




(54) 

(55) 

(56) 

(57) 
(58) 



where m = The foUowing equations can be used for reconstruction of the motion 

of the second body in the inertial frame: 




1 ^ dU 
m2 dX 

X2 = V2, 



(59) 

(60) 
(61) 



i?2 = R2S{n2). 



These equations are equivalent to those given in [12]. However, (59) is not given in [12]. 
Equations (54) though (61) give a complete set of equations for the reduced dynamics and 
reconstruction. Furthermore, they are derived systematically in the context of geometric 
mechanics using proper variational formulas given in (42) through (47). This result can be 
readily generalized for n bodies. 

3.2.2 Equations of motion: Hamiltonian form 

Define the linear momenta r,72 G M^, and the angular momenta 11,112 € as 



Then, the equations of motion can be rewritten in terms of these momenta variables. The 
continuous equations of relative motion for the full two body problem, in Hamiltonian form, 
can be written as 



r = mV, 



72 = mv2, 

n = JrQ = RJiUi, 

U2 = J2^2- 




(62) 
(63) 



n + ri2 X n = -m. 



ri2 + ^^2 X Ha = X X — + M, 



(64) 



X + Q.2XX = — 



m ' 



(65) 
(66) 



R = sin)R - s{n2)R, 
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where m = ^^^^ . The following equations can be used to reconstruct the motion of the 
second body in the inertial frame: 

72 = R2^, (67) 
X2 = —, (68) 

1712 

i?2 = R2S{n2). (69) 



4 Variational integrators 



A variational integrator discretizes Hamilton's principle rather than the continuous equa- 
tions of motion. Taking variations of the discretization of the action integral leads to 
the discrete Euler-Lagrange or discrete Hamilton's equations. The discrete Euler-Lagrange 
equations can be interpreted as a discrete Lagrangian map that updates the variables in 
the configuration space, which are the positions and the attitudes of the bodies. A discrete 
Legendre transformation relates the configuration variables with the linear and angular mo- 
menta variables, and yields a discrete Hamiltonian map, which is equivalent to the discrete 
Lagrangian map. 

In this section, we derive both a Lagrangian and Hamiltonian form of variational integrators 
for the full body problem in inertial and relative coordinates. The second level subscript 
k denotes the value of variables dX t = kh + to for an integration step size /i G R and an 
integer k. The integer A'^ satisfies tf = kN + Iq, so A^ is the number of time-steps of length 
h to go from the initial time to to the final time tf. 



4-1 Inertial coordinates 



Discrete Lagrangian: In continuous time, the structure of the kinematics equations (28), 
(58) and (61) ensure that Ri, R and R2 evolve on SO (3) automatically. Here, we introduce 
a new variable Fi^, G SO (3) defined such that = Ri^Fi^^, i.e. 

Fik = RikRik+1- (70) 

Thus, Fii_ represents the relative attitude between two integration steps, and by requiring 
that Fjj. G S0(3), we guarantee that Ri^ evolves in S0(3). 

Using the kinematic equation Ri = RiS{'D,i), the skew-symmetric matrix S{^i^) can be 
approximated as 

S{n,,) = RlK - Rl^^^'^^^ = liF^k - ^3X3). (71) 

The velocity, be approximated simply by (a^j^^j — a;^^. ) /h. Using these approxima- 

tions of the angular and linear velocity, the kinetic energy of the ith body given in (12) 



18 



can be approximated as 



^11 _ A_ 

1 II Il2 , 1 



where (5) is used. A discrete Lagrangian Lrf(x^,x^^^,Rj^,,F^) is constructed such that it 
approximates a segment of the action integral (21), 

h f 1, X 1. 



i=l 



-^C/(x„RJ-^C/(x,^,,R,^J, (72) 



where x, G (R''^)", R^ G S0(3)", and G (M^)™^ and I G (M^xs^n ^re defined as 
^fc = (a;ife)a;2fe, • • • ,a;nj, Rfc = (^i^ , -^2^ , • • • j-^nj, = (^i^, -^2^, • • • ^-^nj, and I = 
(-^3x3, -^3x3, • ■ ■ , 4x3), respectively. 

This discrete Lagrangian is sclf-adjoint [4], and sclf-adjoint numerical integration methods 
have even order, so we are guaranteed that the resulting integration method is at least 
second-order accurate. 

Variations of discrete variables: The variations of the discrete variables are chosen to 
respect the geometry of the configuration space SE(3). The variation of Xi^ is given by 

where Sxi^, G and vanishes at A; = and k = N. The variation of Ri^ is given by 

SRi^ = Ri^iji^, (73) 

where rji^ G 50 (3) is a variation represented by a skew-symmetric matrix and vanishes at 
A; = and k = N. The variation of Fi^, can be computed from the definition Fj^, = R^^Rj,^^^ 
to give 

5Fi^ = SRj^R'ik+i + RjjRik+i^ 

= -riik^ik + ^ik^ik+i- (74) 
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4.1.1 Discrete equations of motion: Lagrangian form 



To obtain the discrete equations of motion in Lagrangian form, we compute the variation 
of the discrete Lagrangian from (19), (73) and (74), to give 

" 1 1 



i=l 



2 I dxii. '''' dx 



OR 



(75) 



where C/^. = J7(x^,Rj,) denotes the value of the potential sX t = kh + Iq. 
Define the action sum as 



N-l 



(76) 



k=0 



The discrete action sum (5(i approximates the action integral (21), because the discrete 
Lagrangian approximates a segment of the action integral. 

Substituting (75) into (76), the variation of the action sum is given by 

JV-l n 



k=0 i=l ^ 



hdu^ 

2 dr- 



+ tr 
+ tr 



dU, 



k+l 



2 



■I'k 



Using the fact that 5xi^, and "qi^, vanish at A; = and k = N,we can reindex the summation, 
which is the discrete analogue of integration by parts, to yield 



k=i i=i ^ ^^'k ^ 



+ tr 



Vik \ ^ i^ik-^di - JdiFi^-i) + ^^fkQ^ 



Hamilton's principle states that 6<5d should be zero for all possible variations ().x.j^ € 
and ?7ij. G so (3) that vanish at the endpoints. Therefore, the expression in the first brace 
should be zero, and since ryj^. is skew-symmetric, the expression in the second brace should 
be symmetric. Thus, we obtain the discrete equations of motion for the full body problem, 
in Lagrangian form, for i G (1, 2, • • • , n) as 



dxi^ ' 



(77) 
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- [f,^^,M - Jd^Fl^, - JdA. + F^Jd?) = hSiM,^^,), (78) 

-^ifc+i = ^kFiki ('''9) 

where Mi^, e is defined in (24) as 

where rij^,Urip eM}^^ are pth row vectors of Ri,, and respectively. Given the initial 

conditions (a^j,,, , i^i^), we can obtain Xi.^ from (77). Then, _Fj„ is computed from (79), 

and Fi-^ can be obtained by solving the implicit equation (78). Finally, Ri^ is found from 
(79). This yields an update map (xi^, RiQ,Xi-^, Ri-^) t—^ {xi-^, Ri-^, Xi^, Ri^), and this process 
can be repeated. 



4.1.2 Discrete equations of motion: Hamiltonian form 

As discussed above, equations (77) through (79) defines a discrete Lagrangian map that 
updates Xi^. and Ri/^. The discrete Legendre transformation given in (3) and (4) relates the 
configuration variables rcj^, iij^ and the corresponding momenta. This induces a discrete 
Hamiltonian map that is equivalent to the discrete Lagrangian map. The discrete Hamil- 
tonian map is particularly convenient if the initial conditions are given in terms of the 
positions and momenta at the initial time, (x^o, Vig, i^ig, Jljg). 

Before deriving the variational integrator in Hamiltonian form, consider the momenta con- 
jugate to Xi and Ri, namely P^. G and PQ^ G M^^^. Prom the definition (1), F^.L is 
obtained by taking the derivative of L, given in (15), with Xj while holding other variables 
fixed. 

<5i;fP^. =F„.L(x,x,R,J2), 

L(x, X -|- e(5xi, R, fi), 

e=0 



de 
d_ 

Te 



Ti{xi + €dxi,Qi), 



t=0 



where 5xj G 



= Sxi (miXi) , 

)" denotes (0, 0, • • • , dxi, • • • ,0), and Tj is given in (12). Then, we obtain 

Py^ = miVi = ji, (81) 

which is equal to the linear momentum of Bi. Similarly, 

ti[S{SnifPn,] =Fs7,L(x,x,R,n), 

= Ti{xi,ni + edQi), 

= ltr[s{6ni)Jd,sinif + sini)Ja,s{5nif] , 



= -tr:[S{6nifS{Ji^}i)] , 
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where (5) and (13) are used. Now, we obtain 



tr 



0. 



Since 5(r2j) is skew-symmetric, the expression in the braces should be symmetric. This 
imphes that 



(82) 



Equations (81) and (82) give expressions for the momenta conjugate to Xi and Ri. Consider 
the discrete Legendre transformations given in (3) and (4). Then, 



) = - 



Ld(x, +e5xj,,x R,,F 



e=0 



-6x1 



1 . _ ^,!^9u. 



h 



2 dx 



where S^x.i^^ G (M^)" denotes (0, 0, • • • , 6x^^ , ■ ■ ■ ,0). Therefore, we have 

1 \ ^ 

^Xi^k-^d{^k^'^k+i^ ^k^ ^ k) ~ ~ 'g^iy^i-k+i ~ ^ik) ~ 



2dx. 



(83) 



(84) 



Prom the discrete Legendre transformation given in (3), Py^^, = —Dxik^d- Using (81) and 
(84), we obtain 



_ 1 _ hdU^ 



2dxi 



(85) 



Using the discrete Legendre transformation given in (4), P^^ ^^^^ — k+i-^di W6 can derive 
the fohowing equation similarly: 



2dx^ 



(86) 



Equations (85) and (86) define the variational integrator in Hamiltonian form for the 
translational motion. Now, consider the rotational motion. We have 



tr[r7i^Dij,^,L/] = tr 



1 h T dU^ 



(87) 



where the right side is obtained by taking the variation of with respect to Ri^, while 
holding other variables fixed. Since r)ii^ is skew-symmetric. 



(88) 



where Mj. G is defined in (80). 
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Prom the discrete Legendre transformation given in (3), Pcii^ = —^Rik^d, we obtain the 
following equation by using (82) and (88), 

Si^iJ = I {F^,Jd. - Jd,Fl) - ^5(M,J. (89) 

Using the discrete Legendre transformation given in (4), Pq. = Dr. ^_^_j^La, we can obtain 
the following equation: 

S{^i,+,) = {Fi.Jd, - Ja^Fl) F,^ + |5(M,,^J. (90) 

By using (10) and substituting (89), we can reduce (90) to the following equation in vector 
form. 

n.,+, = F^n,, + |i^^M,, + |m,,^, . (91) 

Equations (89) and (91) define the variational integrator in Hamiltonian form for the 
rotational motion. 

In summary, using (85), (86), (89) and (91), the discrete equations of motion for the full 
body problem, in Hamiltonian form, can be written for i G (1, 2, • ■ ■ , n) as 

h _ h'^ dU^ 



h 




mi 




hdU^ 


h dU,^, 


2dxi^ 


2 8t 



h 



hS{Ui, + --^Mi,) = Fija, - Jd,Fl, (94) 



2 



„'r , , h 



n,,+, = Flli,^ + -F^M,^ + '-M,^^^., (95) 

Rik+i = RikFik- (96) 

Given (xjg, 7jg, nj^), we can find Xj^from (92). Solving the implicit equation (94) yields 
Fjo, and Ri^ is computed from (96). Then, (93) and (95) gives 7^^, and Ilj^. This defines 
the discrete Hamiltonian map, (xjg , , Ri^ , lii^ ) ^ {xi^ , 7^^ , Ri^ , 11,^ ) , and this process can 
be repeated. 

4-2 Relative coordinates 

In this section, we derive the variational integrator for the full two body problem in relative 
coordinates by following the procedure given before. This result can be readily generalized 
to n bodies. 

Reduction of discrete variables: The discrete reduced variables are defined in the same way 
as the continuous reduced variables, which are given in (34) through (39). We introduce 
G S0(3) such that R,^^ = R^.^Ri,^, = F^F.R^. i.e. 

F,=R,F^,Rl. (97) 
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Discrete reduced Lagrangian: The discrete reduced Lagrangian is obtained by expressing 
the original discrete Lagrangian given in (72) in terms of the discrete reduced variables. 

Prom the definition of the discrete reduced variables given in (34) and (38) , we have 

= R2, {F2,(X,+, + ^2,+ J - {X, +X2J} , (98) 
X2,^, X2, = R2, {F2,X2,^-, -X2,]. (99) 

From (71), S{Q,i^) and S{Q,2k) are expressed as 

'5(^^iJ = ;^(i^i.-/3x3), 

= \Rl{Fu-h>.^)R,, (100) 
'5(^^2j = ;^(F2,-/3x3). (101) 
Substituting (98) through (101) into (72), we obtain the discrete reduced Lagrangian. 

ku = h{X^ , X^_^^ , X2^ , X2,^^^ ,Rk^Fk^ ) 

= ^mi \\F2,{X,^,+X2,^,) - {X,+X2,)f + ^m2 \\F2,X2,^, - 

+ ^tr[(/3x3 - F,)JaR,] + ^tr[(/3x3 - F2M - ^U{X^,R,) - ^U{X^^,, R,^,), 

(102) 

where JdR^ ^ K^^^ is defined to be JdR^ = Rk'^diR^-: which gives the nonstandard moment 
of inertia matrix of the first body with respect to the second body fixed frame at t = kh+to- 

Variations of discrete reduced variables: The variations of the discrete reduced variables 
can be derived from those of the original variables. The variations of R^, X^,, and are 
the same as given in (42), (43), and (74), respectively. The variation of is computed in 
A.2. 

In summary, the variations of discrete reduced variables are given by 



SR,=riA-mA, (103) 

6X^=X,-m,X„ (104) 

SF, = -m,F, + F2,rj,^,FlF, + F, {-rj, + m,) , (105) 

^^2, =X2,-r?2,X2„ (106) 

SF2, = -V2A+F2,m,+r- (107) 



These Lie group variations are the main elements required to derive the variational inte- 
grator equations. 

4.2.1 Discrete equations of motion: Lagrangian form 

As before, we can obtain the discrete equations of motion in Lagrangian form by computing 
the variation of the discrete reduced Lagrangian which, by using (103) through (107), is 
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given by 



1 



h 

+ [mi{X^ + X2,) - miF2,(X,+i + + m2X2, - ^2^2^X2,^^ 



h J, dUf. h 



h 



dX,_ 



dX„ 



h 

77X 



fe+i 



T- 



h 



fc+1 fc+1 



(108) 



The action sum expressed in terms of the discrete reduced Lagrangian has the form 



N-l 



6. 



(109) 



fe=0 



The discrete action sum 0^ approximates the action integral (49) , because the discrete La- 
grangian approximates a piece of the integral. Using the fact that the variations Xk 1 X2u 1 Vk ' V2k 
vanish at A; = and k = N, the variation of the discrete action sum can be expressed as 



7V-1 , 

= J2 TT^r - miFl_^{X,_, + X2,_J + 2mi(X, + X2J 
k=l ^ 

-miF2,(X,^,+X2,^J-/i2|^| 

N-l ^ . 
k=l ^ 

- m2Fl^_^X2,_, + 2m2X2, - 7712^2, X2,^, | 

+ E il [-Pl-.P.-.K-.Jd,Rl_F2,_,+F,RJa,Rl)-hR,^ 
k=i ^ ^ 

N—l |- f ^ r^j-j. J- j^^^ J-,, -| 

+ ^ tr 



k=l 



mk 1 1 {-JdA_, + F2jd,) + ^^^1^ + ^^^1^ 



Prom Hamilton's principle, (5(5^ should be zero for all possible variations x. 1 X2k ^ ^ 
and ?7fc,?72j. G so (3) which vanish at the endpoints. Therefore, in the above equation, the 

expressions in the first two braces should be zero, and the expressions in the last two braces 
should be symmetric since r/^ , 772. are skew-symmetric. After some simplification, we obtain 
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the discrete equations of relative motion for the full two body problem, in Lagrangian form, 
as 

F,^X,^, - 2X, + Fl_^X,_, = -^1^, (110) 
Fk+iJdRk+i - JdRk+iFk+i = ^^k i^kJdRk - JdRk^k) ^^k " h'^S{Mk+i), (111) 

^ ^^yi / \ 2 ____ ^^Cy 2 / X 

F2k+iJd2 - -^^2-^2^+1 = -^2fc (-^2fc-'d2 - JdiF^k) F% + ^ ^k+l ^ 1- ^ 'S'(Mfc+i), 

(112) 

Rk+.=FlFA, (113) 

R2,+,=R2A. (114) 

It is natural to express equations of motion for the second body in the inertial frame. 

^^k+i - 2^2, + X2,_, = ^^fc (115) 

Given {X^, R^^, R2q, X-^, R^, R2j^), we can determine and from (113) and (114). Solv- 
ing the implicit equations (111) and (112) gives and F2^. Then X^, R^ and R22 are 
found from (110), (113) and (114), respectively. This yields the discrete Lagrangian map 
(Xg , i?g , i?2o ) ^1 , -^1 , -^21 ) I— i^i, Rn R21, X^, R^, R22) and this process can be repeated. 
We can separately reconstruct X2^. using (115). 



4-2.2 Discrete equations of motion: Hamiltonian form 

Using the discrete Legendre transformation, we can obtain the Hamiltonian map, in terms 
of reduced variables, that is equivalent to the Lagrangian map given in (110) through 
(115). We will only sketch the procedure as it is analogous to the approach of the previous 
section. First, we find expressions for the conjugate momenta variables corresponding to 
(81) and (82). We compute the discrete Legendre transformation by taking the variation of 
the discrete reduced Lagrangian as in (83) and (87). Then, we obtain the discrete equations 
of motion in Hamiltonian form using (3) and (4). 

The discrete equations of relative motion for the full two body problem, in Hamiltonian 
form, can be written as 

X^, = Fl (x^ + h^- , (116) 

'=+1 k m 2m dX^ J' ^ ' 

^--^^^i^^-2 9x:j-2ax:;:' ^'''^ 



n,+i = Fl ( H, - - ^M,^,, (118) 



n2.,. ^ Fl (h. H- X ^ + ^M,) + ^X.,, X ^ + ^M.,,, (119) 

Rk^.=FlF,R,, (120) 



hS[Iik- J = FJaR, - Jm,F^, (121) 
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hS (U2, + X ^ + 1^.) = F2,Ja, - Jd^Fl- (122) 

It is natural to express equations of motion for the second body in the inertial frame for 
reconstruction: 

X2 = X2, + + -^R,^, (123) 



h da h dU^,^ , ^ 

72,+, = 72, + 2^'=a^ + 2^^+i aX^' ^^^^^ 



i?2,+i=i?2,i^2,. (125) 

Given {R^ , , , Fg , n2(j ) , we can determine and by solving the implicit equations 
(121) and (122). Then, and arc found from (116) and (120), respectively. After that, 
we can compute F^, U^, and 112^ from (117), (118) and (119). This yields a discrete Hamil- 
tonian map {R^ , X^ , IIq , Fg , n2o ) 1-^ {R^ ) ^1 ) ) ) ) , and this process can be repeated. 
X2ki 72, and i?2, can be updated separately using (123), (124) and (125), respectively, for 
reconstruction. 



If..3 Numerical considerations 



Properties of tfie variational integrators: Variational integrators exhibit a discrete analogue 
of Noether's theorem [17], and symmetries of the discrete Lagrangian result in conservation 
of the corresponding momentum maps. Our choice of discrete Lagrangian is such that it in- 
herits the symmetries of the continuous Lagrangian. Therefore, all the conserved momenta 
in the continuous dynamics are preserved by the discrete dynamics. 

The proposed variational integrators are expressed in terms of Lie group computations [5] . 
During each integration step, Fi^ G SO (3) is obtained by solving an implicit equation, and 
Ri^ is updated by multiplication with Fi^,. Since SO (3) is closed under matrix multiplica- 
tion, the attitude matrix Ri^^^ remains in S0(3). We make this more explicit in section 
4.4 by expressing Fi^_ as the exponential function of an element of the Lie algebra so (3). 

An adjoint integration method is the inverse map of the original method with reversed 
time-step. An integration method is called self-adjoint or symmetric if it is identical with 
its adjoint; a self-adjoint method always has even order. Our discrete Lagrangian is chosen 
to be self-adjoint, and therefore the corresponding variational integrators are second-order 
accurate. 

Higher-order methods: While the numerical methods we present in this paper are second- 
order, it is possible to apply the symmetric composition methods, introduced in [26], to 
construct higher-order versions of the Lie group variational integrators introduced here. 
Given a basic numerical method represented by the flow map the composition method 
is obtained by applying the basic method using different step sizes, 

= ^Xsh o • • • o $Ai/i, 
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where Ai,A2,-- - ,As G M. In particular, the Yoshida symmetric composition method for 
composing a symmetric method of order 2 into a symmetric method of order 4 is obtained 
when s = 3, and 

1 2V3 



Alternatively, by adopting the formalism of higher-order Lie group variational integrators 
introduced in [11] in conjunction with the Rodrigues formula, one can directly construct 
higher-order generalizations of the Lie group methods presented here. 

ReducUon of orthogonality loss due to roundoff error: In the Lie group variational inte- 
grators, the numerical solution is made to automatically remain on the rotation group by 
requiring that the numerical solution is updated by matrix multiplication with the expo- 
nential of a skew symmetric matrix. 

Since the exponential of the skew symmetric matrix is orthogonal to machine precision, 
the numerical solution will only deviate from orthogonality due to the accumulation of 
roundoff error in the matrix multiplication, and this orthogonality loss grows linearly with 
the number of timesteps taken. 

One possible method of addressing this issue is to use the Baker-Campbell-Hausdorff (BCH) 
formula to track the updates purely at the level of skew symmetric matrices (the Lie 
algebra). This allows us to find a matrix C{t), such that, 

ex.p{tA) exp{tB) = expC(t). 

This matrix C{t) satisfies the following differential equation, 

C = A + B + ^[A-B,C]+^^ad'l,{A + B), 

k>2 

with initial value C(0) = 0, and where B^ denotes the Bernoulli numbers, and adc^ = 
[C, A]=CA- AC. 

The problem with this approach is that the matrix C{t) is not readily computable for 
arbitrary A and B, and in practice, the series is truncated, and the differential equation is 
solved numerically. 

An error is introduced in truncating the series, and numerical errors arc introduced in 
numerically integrating the differential equations. Consequently, while the BCH formula 
could be used solely at the reconstruction stage to ensure that the numerical attitude always 
remains in the rotation group to machine precision, the truncation error would destroy the 
symplecticity and momentum preserving properties of the numerical scheme. 

However, by combining the BCH formula with the Rodrigues formula in constructing the 
discrete variational principle, it might be possible to construct a Lie group variational 

integrator that tracks the reconstructed trajectory on the rotation group at the level of a 
curve in the Lie algebra, while retaining its structure-preservation properties. 
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4-4 Computational approach 



The structure of the discrete equations of motion given in (78), (94), (111), (112), (121), 
and (122) suggests a specific computational approach. For a given g G M^, we have to solve 
the following Lyapunov-like equation to find F G SO (3) at each integration step. 

FJd - JdF'^ = S{g). (126) 

We now introduce an iterative approach to solve (126) numerically. An clement of a Lie 
group can be expressed as the exponential of an clement of its Lie algebra, so F G SO (3) 
can be expressed as an exponential of S{f) € so (3) for some vector / G R^. The exponential 
can be written in closed form, using Rodrigues' formula, 

F = e^^f\ 



/3x3 + ^5(/) + l^^W. (127) 



Substituting (127) into (126), we obtain 



where (9) and (13) arc used. Thus, (126) is converted into the equivalent vector equation 
g = G{f), where G : is 

Gif) = 'Wjf+'-^^f><Jf- 



We use the Newton method to solve g = G{f), which gives the iteration 

fi+i = fi + VGifi)-^ {g - G{fi)) . (128) 

We iterate until — G(/i)|| < e for a small tolerance e > 0. The Jacobian VG(/) in (128) 
can be expressed as 



^G{f) = 11-^11 IIJII -siJ^ 11/ II jjjT ^ sin||, 



^ sin 11/11 11/11 -2(1- cos ^ 



+ ^-^^{-s{Jf) + s{f)J}. 



Numerical simulations show that 3 or 4 iterations are sufficient to achieve a tolerance of 
e= 10-^^ 



5 Numerical simulations 



The variational integrator in Hamiltonian form given in (116) through (125) is used to 
simulate the dynamics of two simple dumbbell bodies acting under their mutual gravity. 
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5.1 Full body problem defined by two dumbbell bodies 



Each dumbbell model consists of two equal rigid spheres and a massless rod as shown in 
Fig. 3. The gravitational potential of the two dumbbell models is given by 



U{X, R) 



2 

V - Gmim2/4: 
' \\X + p2+Rpi^ 

p,q=l II ' P ' 1 



where G is the universal gravitational constant, € M is the total mass of the ith dumb- 
bell, and Pip G M'^ is a vector from the origin of the body fixed frame to the pth sphere of 
the ith dumbbell in the ith body fixed frame. The vectors pi-^ = [/j/2, 0, 0]-'", 
where is the length between the two spheres. 



Pi2 



-Pii, 





Fig. 3. Dumbbell model of the full body problem 

Normalization: Mass, length and time dimensions are normalized as follows, 

mi 

mi = — , 
m 

X 



G{mi + 7712) 

where m = JUiIB^ and / is chosen as the initial horizontal distance between the center of 
mass of the two dumbbells. The time is normalized so that the orbital period is of order 
unity. Over-bars denote normalized variables. We can expresses the equations of motion in 
terms of the normalized variables. For example, (54) can be written as 

dX 

where ' denotes a derivative with respect to i. The normalized gravitational potential and 
its partial derivatives are given by 
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Conserved quantities: The total energy E is conserved: 
E=^mi \\V + V2f + ^m2\\V2f 

+ ^tr[sin)Ja^s{nf] + ^tr[s{n2)Jd,s{n2f] + u{x,r). 

The total linear momentum 7^ € M"^, and the total angular momentum about the mass 
center of the system ttt S M^, in the inertial frame, are also conserved: 

7T = R2 {mi {V + V2) + m2V2} , 

TTT = R2 {mX XV + JrQ. + J2^2} ■ 



5.2 Simulation results 



The properties of the two dumbbell bodies are chosen to be 

mi = 1.5, h = 0.25, Ji = diag 0.0004 0.0238 0.0238 

7712 = 3, I2 = 0.5, J2 = diag 0.0030 0.1905 0.1905 

The mass and length of the second dumbbell are twice that of the first dumbbell. The 
initial conditions are chosen such that the total linear momentum in the inertial frame is 
zero and the total energy is positive. 



^0 = 



1 0.3 
-0.33 -0.1 



V2o 



1 
-0.33 



^2o = 



9 




, = 73x3, 

, i?2o = -^3x3- 



Simulation results obtained using the Lie group variational integrator are given in Fig. 4 
and Fig. 5. Fig. 4 shows the trajectory of the two dumbbells in the inertial frame. Fig. 5(a) 
shows the evolution of the normalized energy, where the upper figure gives the history of the 
translational kinetic energy and the rotational kinetic energy, and the lower figure shows 
the interchange between the total kinetic energy and the gravitational potential energy. 
Fig. 5(b) shows the evolution of the theoretically conserved quantities, where the upper 
figure is the history of the total energy, and the lower figure is the error in the rotation 
matrix. 

Initially, the first dumbbell rotates around the vertical 63 axis, and the second dumb- 
bell does not rotate. Since the angular velocity of the first dumbbell is relatively large. 
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Fig. 4. Trajectory in the inertial frame 




(a) Interchange of energy (b) Conserved quantities 

Fig. 5. Lie group variational integrator 

the rotational kinetic energy initially exceeds the translational kinetic energy. As the two 
dumbbells orbit around each other, the second dumbbell starts to rotate, the rotational 
kinetic energy increases, and the translational kinetic energy decreases slightly for about 
6 normalized units of time. At 9 units of time, the distance between the two dumbbells 
reaches its minimal separation, and the potential energy is transformed into kinetic energy, 
especially translational kinetic energy. After that, two dumbbells continue to move apart, 
and the translational energy and the rotational energy equalize. (A simple animation of 
this motion can be found at ^http : //www . umich . edu/~ tylee.) This shows some of the 
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interesting dynamics that the fuh body problem can exhibit. The non-trivial interchange 
between rotational kinetic energy, translational kinetic energy, and potential energy may 
yield complicated motions that cannot be observed in the classical two body problem. 

The Lie group variational integrator preserves the total energy and the geometry of the 
configuration space. The maximum deviation of the total energy is 2.6966 x 10""^, and the 
maximum value of the rotation matrix error / — R'^R\\ is 2.8657 x 10^^^. 




5 10 15 5 10 15 



t t 

(a) Interchange of energy (b) Conserved quantities 

Fig. 6. Runge-Kutta method 

As a comparison. Fig. 6 shows simulation results obtained by numerically integrating the 
continuous equations of motion (62)- (69) using a standard Runge-Kutta method. The ro- 
tational and the translational kinetic energy responses are similar to those given in Fig. 
5 prior to the close encounter. However, it fails to simulate the rapid interchange of the 
energy near the minimal separation of the two dumbbells. The deviation of the total energy 
is relatively large, with a maximum deviation of 1.1246 x 10~^. Also, the energy transfer 
is quite different from that given in Fig. 5(a). The Runge-Kutta method does not preserve 
the geometry of the configuration space, as the discrete trajectory rapidly drifts off the 
rotation group to give a maximum rotation matrix error of 2.2435 x 10~^. As the gravity 
and momentum between the two dumbbells depend on the relative attitude, the errors in 
the rotation matrix limits the applicability of standard techniques to long time simulations. 



6 Conclusions 

Eight different forms of the equations of motion for the full body problem are derived. 
The continuous equations of motion and variational integrators are derived both in inertial 
coordinates and in relative coordinates, and each set of equations of motion is expressed 
in both Lagrangian and Hamiltonian form. The relationships between these equations of 
motion are summarized in Fig. 7. This commutative cube was originally given in [6] . In the 
figure, dashed arrows represent discretization from the continuous systems on the left face of 
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the cube to the discrete systems on the right face. Vertical arrows represent reduction from 
the full (incrtial) equations on the top face to the reduced (relative) equations on the bottom 
face. Front and back faces represent Lagrangian and Hamiltonian forms, respectively. The 
corresponding equation numbers are also indicated in parentheses. 

Hamilton Discrete Hamilton 

(29)-(32) ^ (92)-(96) 




Reduced EL Reduced DEL 

(54)-(61) Di^etizaiicni^ ^ (110)-(115) 

Fig. 7. Commutative cube of the equations of motion 

It is shown that the equations of motion for the full body problem can be derived system- 
atically, using proper Lie group variations, from Hamilton's principle. The proposed varia- 
tional integrators preserve the momenta and symplcctic form of the continuous dynamics, 
exhibit good energy properties, and they also conserve the geometry of the configuration 
space since they are based on Lie group computations. The main contribution of this paper 
is the combination of variational integrators and Lie group computations, developed for the 
full body problem. Hence, the resulting numerical integrators conserve the first integrals 
as well as the geometry of the configuration space of the full body dynamics. 

A Appendix 

A.l Variations of reduced variables 

The variations of the reduced variables given in (43) through (47) are derived in this section. 
The variations of the reduced variables can be obtained from the definitions of the reduced 
variables, and the variations of the original variables. 

The variation of X = i?J(xi — X2) is given by 

5X = 5R^{x\ — X2) + R2{6xi — 6x2)- 
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s{dn) = ^ 



e=0 



Substituting (16) into the above equation, we obtain 

SX = -rj2R2ixi - X2) + R2{Sxi - 6x2), 
= -ri2X + X, 

where the reduced variation x '■ [^O; '—>■ is defined to be x = R2{^xi — 6x2)- 
From the definition of = RQi and (10), S{SQ,) is given by 

R'S{n\)R'^, 

) 

= 6RS{ni)R^ + RS{6ni)R^ + RS{ni)6R^. 
Substituting (42) and (17) into the above equation, we obtain 

s{6n) = {r]- m} RS{ni)R^ + R{m + s{ni)r]i - mSi^i)} R^ 

+ RS{ni)R^{-r] + r]2}, 
= {t1- m} SiRQi) + RmR'^ + S{Rni)RriiR^ - Rr]iR^S{Rni) 
+ SiRni){-r] + r]2}. 

Since rj = RijiR^ and = RQi, the above equation reduces to 

S{6n) = -mSin) + RrjiR^ + S{n)rj2- (A.l) 

Prom the definition of i? = R2R1, R is given by 

R = R^ Ri + R^Ri, 

= -S{n2)R + S{n)R. (A.2) 

Then, rj can be written as 

■f] = RriiR^ + RriiR^ + Rr]iR^ , 
= Rr)iR^ + {S{n) - S{n2)}r] - i] {S{n) - ^(^2)} ■ (A.3) 

Substituting (A.3) into (A.l), we obtain S{5Vl) in terms of ?7,r/2 as 

S{5Q.) = 77 - S{^)ri + riS{a) + S{a)r]2 - mS{^) + S{Vt2)r] - r?5(Q2), 

which is equivalent to (44). 

The variation oiV = R^{x\ — X2) is given by 

5V = 5R^{xi - ±2) + R2{5xi - 8x2), 

= -mV + R^iSxi -6x2). (A.4) 

Prom the definition of % = i?J((5xi — 6x2), X is given by 

X = R^iSxi - 6x2) + R^iSxi - 6x2), 
= ~S{n2)x + Rl{5xi - 8x2). (A.5) 
Substituting (A.5) into (A.4), we obtain 

5v = -mv + X + Sin2)x, 

which is equivalent to (45). The variation SV2 can be derived in the same way, and 5(5^2) 
is given in (17). 
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A. 2 Variations of discrete reduced variables 

The variation of the reduced variables 6F^ given in (105) is derived in this section. Prom 
(74) and (97), the variation 6Fi^, is written as 

where r?, G 5o(3) is defined as r,, = R,m,Rl- Since F^R^Rl^^ = F^R^iR^ F^ F^,) = F2,, 
we have 

•^^i. = (-^.^/c + F2,V,+,FlF^) R^. 
Then, the variation SF^, is given by 

SF^ = SR,F,^R^ + RJF.^R^ + R.F.^SR^, 

which is equivalent to (105). 
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